module func_drv_test; 
use numeric_deriv;
use numeric_func;
implicit none;

contains 

logical function func_drv_test_all() result(ok);  implicit none; 
  call dbg_info("-->func_drv_test_all()"); 
  ok = .TRUE.
  ok = new_func_drv_test();
  !ok = vec_dot_test_speed();
  return
end function

logical function new_func_drv_test() result(ok);  implicit none; 
  integer :: nv;
  real(R14) :: minv, maxv, res, MAX_DRV_ERR = 1.d-18;
  real(R14), dimension(200) :: px, py, py2, pdy, pf;      
  real(R14), dimension(3) :: coeff, dc;      
  real(R14), dimension(6) :: c6, dc6;      
  real(R14), dimension(2) :: c2, dc2;      
  call dbg_info("-->new_func_drv_test()"); 
  ok = .TRUE.

  nv = 9;  minv = 0.; maxv = 1.;
!  px => new_step_grid(minv, maxv, nv);
  call load_step_grid(px, 1, nv, minv, maxv);
  call assert_max_zero("9.px(1)=minv", px(1) - minv);
  call assert_max_zero("9.px(2)=minv", px(2) - (maxv+minv)/8.);
  call assert_max_zero("9.px(nv)=maxv", px(nv) - maxv);
  
  coeff = (/1, 0, 0/);  ! y=1
  call load_func_polynom(py, px, coeff);  
  call load_func_deriv(pdy, px, py);
  dc = (/0, 0, 0/);  ! y=0
!  pf => new_vec_polynom(px, dc);  
  call load_func_polynom(pf, px, dc);  
  call assert_max_diff_vec("pdy(1)", pdy, pf, MAX_DRV_ERR);

  coeff = (/D1, D1, D0/);  ! y=1+x
  py => new_vec_polynom(px, coeff);  
  pdy => new_func_drv(px, py);
  res = D1;
  call assert_max_diff("pdy(1)", pdy(1), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(2)", pdy(2), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(3)", pdy(3), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(4)", pdy(4), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(5)", pdy(5), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(6)", pdy(6), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(7)", pdy(7), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(8)", pdy(8), res, MAX_DRV_ERR);
  call assert_max_diff("pdy(9)", pdy(9), res, MAX_DRV_ERR);
  call del_vec(py);
  call del_vec(pdy);

  coeff = (/D0, D0, D05/);  
  py => new_vec_polynom(px, coeff);  
  pdy => new_func_drv(px, py);
  dc = (/D0, D1, D0/);  
  pf => new_vec_polynom(px, dc);  
  call assert_max_diff_vec("y=0.5*x^2", pdy, pf, MAX_DRV_ERR);
  call del_vec(py);
  call del_vec(pdy);
  call del_vec(pf);

  coeff = (/D0, D0, -D05/);  
  py => new_vec_polynom(px, coeff);  
  pdy => new_func_drv(px, py);
  dc = (/D0, -D1, D0/); 
  pf => new_vec_polynom(px, dc);  
  call assert_max_diff_vec("y=-0.5*x^2", pdy, pf, MAX_DRV_ERR);
  call del_vec(py);
  call del_vec(pdy);
  call del_vec(pf);

  c6 = (/-D01, -D05, -D05, D0, D0, D0/);  
  py => new_vec_polynom(px, c6);  
  pdy => new_func_drv(px, py);
  dc6 = (/-D05, -D1, D0, D0, D0, D0/);  
  pf => new_vec_polynom(px, dc6);  
  call assert_max_diff_vec("y=-0.1-0.5*x-0.5*x^2", pdy, pf, ERR14);
  call del_vec(py);
  call del_vec(pdy);
  call del_vec(pf);

  c6 = (/-D02, D05, -D05, D2, -1D0/4D0, D0/);  
  py => new_vec_polynom(px, c6);  
  pdy => new_func_drv(px, py);
  dc6 = (/D05, -D1, D6, -D1, D0, D0/);  
  pf => new_vec_polynom(px, dc6);  
  call assert_max_diff_vec("y=...x^4/4", pdy, pf, ERR14);
  call del_vec(py);
  call del_vec(pdy);
  call del_vec(pf);

  py => new_vec_polynom(px, (/D0, D2/));  ! y=2x
  py2 => new_func_exp(py);  ! y=exp(2x)
  pdy => new_func_drv(px, py2);  ! y'=[exp(2x)]'=2exp(2x)
  pf => new_vec_from(py2);  
  pf = py2 * 2;  
  call assert_max_diff_vec("y=exp(2x)", pdy, pf, 2 * ERR5);
  call del_vec(py);
  call del_vec(py2);
  call del_vec(pdy);
  call del_vec(pf);

  ! NEW px
  nv = 101;  minv = -D1; maxv = D1;
  call del_vec(px);
!  px => new_step_grid(minv, maxv, vsize);
  px => new_vec(nv);  
  call load_step_grid(px, 1, nv, minv, maxv);
  py => new_vec_polynom(px, (/D0, D2/));  ! y=2x
  py2 => new_func_exp(py);  ! y=exp(2x)
  pdy => new_func_drv(px, py2);  ! y'=[exp(2x)]'=2exp(2x)
  pf => new_vec_from(py2);  
  pf = py2 * 2;  
  call assert_max_diff_vec("y=exp(2x)", pdy, pf, ERR11);
  call del_vec(py);
  call del_vec(py2);
  call del_vec(pdy);
  call del_vec(pf);

  c6 = (/-D02, D05, -D05, D2, -1D0/4D0, D0/);  
  py => new_vec_polynom(px, c6);  
  pdy => new_func_drv(px, py);
  dc6 = (/D05, -D1, D6, -D1, D0, D0/);  
  pf => new_vec_polynom(px, dc6);  
  call assert_max_diff_vec("y=...x^4/4", pdy, pf, ERR12);
  call del_vec(py);
  call del_vec(pdy);
  call del_vec(pf);


  call del_vec(px);
  return
end function

end module 